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Abstract 

A formulation for the discontinuous Galerkin (DG) method that leads to solutions using the 
differential form of the equation (as opposed to the standard integral form) is presented. The formulation 
includes (a) a derivative calculation that involves only data within each cell with no data interaction 
among cells, and (b) for each cell, corrections to this derivative that deal with the jumps in fluxes at the 
cell boundaries and allow data across cells to interact. The derivative with no interaction is obtained by a 
projection, but for nodal-type methods, evaluating this derivative by interpolation at the nodal points is 
more economical. The corrections are derived using the approximate (Dirac) delta functions. The 
formulation results in a family of schemes: different approximate delta functions give rise to different 
methods. It is shown that the current formulation is essentially equivalent to the flux reconstruction (FR) 
formulation. Due to the use of approximate delta functions, an energy stability proof simpler than that of 
Vincent, Castonguay, and Jameson (2011) for a family of schemes is derived. Accuracy and stability of 
resulting schemes are discussed via Fourier analyses. Similar to FR, the current formulation provides a 
unifying framework for high-order methods by recovering the DG, spectral difference (SD), and spectral 
volume (SV) schemes. It also yields stable, accurate, and economical methods. 

1. Introduction 

For conservation laws, e.g., those in the field of Computational Fluid Dynamics (CFD), low-order 
methods are generally robust and reliable; as a result, they are routinely employed in practical 
calculations. For the same computing cost, high-order methods can provide considerably more accurate 
solutions, but they are more complicated and less robust. The need to improve and develop new high- 
order methods with favorable properties has attracted the interest of many researchers. 

The discontinuous Galerkin (DG) method is currently among the most widely used high-order 
numerical methods for solving conservation laws on unstructured meshes. It was introduced for the 
neutron transport equation by Reed and Hill (1973), analyzed by LaSaint and Raviart (1974) and 
developed and made popular for fluid dynamics equations by Cockbum, Shu, Bassi, Rebay, and others 
(Cockbum, Kamiadakis, and Shu 2000, Bassi and Rebay 1997a,b, 2000, Cockbum and Shu 2005, 2009, 
Shu 2012, and the references therein). Efficient DG schemes storing data at nodal points kn own as nodal 
DG methods were elaborated by Hesthaven and Warburton (2008). 

Alternative approaches to high-order accuracy employing the differential form (as opposed to DG 
which employs the integral form) have been proposed. Kopriva and Kolias (1996) pioneered this 
approach with the staggered-grid spectral method on quadrilateral meshes. The approach was extended to 
triangular meshes by Liu, Vinokur, and Wang (2006) and named spectral difference (SD), and solutions 
for a wide range of problems were obtained by Wang et al. (2007) and Liang et al. (2009a and 2009b). 
Another class of schemes called spectral volume (SV) based on the idea of subdividing each cell into 
subcells or control volumes in a structured manner was proposed by Wang, Zhang, and Liu (2004). 
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Recently, an approach to high-order accuracy with the advantage of simplicity and economy called 
flux reconstruction (FR) was introduced by the author (Huynh 2007, 2009). The approach amounts to 
evaluating the derivative of a discontinuous piecewise polynomial function by employing its 
straightforward derivative estimate together with a correction that accounts for the jumps at the interfaces. 
The FR framework unifies several existing schemes: with appropriate choices of correction terms, it 
recovers DG, SD, as well as SV, and the FR versions are generally simpler and more economical than the 
original versions. In addition, the approach results in numerous new methods that are stable and super 
accurate, i.e., more accurate than expected (also known as super convergent). Extensions of FR to 
unstructured triangular meshes was provided by Wang and Gao (2009), to the 2D Navier-Stokes 
equations on meshes of mixed elements by Gao and Wang (2009) and Gao et al. (2013), to the 3D Euler 
and Navier-Stokes equations on mixed meshes by Haga et al. (2010, 201 1) and Wang et al. (2011), and to 
dynamic meshes by Yu et al. (2012). Whereas FR approach is also known as CPR (Correction Procedure 
via Reconstruction), due to the usage in references closely related to this work, the name FR is employed. 

A mathematical foundation for the FR approach was recently provided by Jameson (2010), who 
proved that a particular SD scheme (recovered via FR) is energy-stable for ID linear advection. Vincent, 
Castonguay, and Jameson (2011a) subsequently extended this result, and proved that a one-parameter 
family of FR methods is energy-stable for linear advection. This family, referred to as Energy Stable Flux 
Reconstruction (ESFR), was extended to linear advection on 2D triangular grids by Castonguay, Vincent 
and Jameson (2012), to linear advection-diffusion in ID by Castonguay et al. (2013), and to linear 
advection-diffusion on 2D triangular grids by Williams et al. (2011) and Williams et al. (2013). Von 
Neumann analyses for these schemes were investigated by Vincent, Castonguay and Jameson (2011b) and 
nonlinear stability was discussed by Jameson, Vincent and Castonguay (2012). The ESFR methods were 
shown to be equivalent a filtered DG scheme by Allaneau and Jameson (2012). 

A review of recent developments for the FR methods was presented by Huynh, Wang, and Vincent 
(2014). An open-source Python based framework for solving advection-diffusion type problems on 
streaming architectures using the FR approach can be found at www.pyff.org/. 

In this paper, a formulation for the discontinuous Galerkin (DG) methods that leads to solutions using 
the differential form of the equation is presented. The formulation includes (a) a derivative calculation 
that involves only data within each cell with no data interaction among cells and (b) for each cell, 
corrections to this derivative that deal with the jumps in fluxes at the cell boundaries and allow data 
across cells to interact. The derivative with no interaction is obtained by a projection, but for nodal-type 
methods, evaluating of this derivative by interpolation at the nodal points is more economical. The 
corrections are derived using the approximate (Dirac) delta functions. The formulation results in a family 
of schemes: different approximate delta functions give rise to different methods. It is shown that the 
current formulation is essentially equivalent to the flux reconstruction (FR) formulation. Due to the use of 
approximate delta functions, an energy stability proof simpler than that of Vincent, Castonguay, and 
Jameson (2011) for a family of schemes is derived. Accuracy and stability of resulting schemes are 
discussed via Fourier analyses. 

The formulations to be discussed are basic in nature. Therefore, in the hope of attracting the interest of 
researchers who are not familiar with (some of) these approaches, this paper is written in an essentially 
self-contained manner. It is organized as follows. The DG method is reviewed in Section 2. Section 3 
introduces new strong forms for DG via the use of the approximate Dirac delta functions. The FR 
schemes, obtained by integrating the strong forms of Section 3, are discussed in Section 4. An energy 
stability proof for a family of schemes is presented in Section 5. Section 6 contains representative FR 
schemes and their properties obtained by Fourier analysis. Conclusions and discussion are presented in 
Section 7. A brief discussion of Fourier stability and accuracy analyses can be found in the Appendix. 
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2. Discontinuous Galerkin Formulations 


2.1. Preliminaries. Consider the conservation law 


u t + fx = 0 ( 2 . 1 ) 

with initial condition u(x, 0) = u init (x). The solution u is assumed to be periodic or of compact support 
so that boundary conditions are trivial. The flux / depends on u. With a(u) = df /du, the above can be 
cast in nonconservation form: 


u t + au x = 0 . 


( 2 . 2 ) 


Let the domain of calculation fi be divided into possibly nonuniform cells or elements £), j = 1, 2 , ... 
Denote the center of Ej by Xj and its width by hj. Instead of dealing with the global elements, it is more 
convenient to deal with the local (or reference) element / = [—1, 1]. With ^ varying on / and x on £), the 
linear function mapping / onto £)■ and its inverse are 

x(<f) = Xj + ^hj/2 and ^(x) = 2(x — Xj)/hj . (2.3) 

A function ry (x) on Ej results in a function on / denoted by ry (^) for simplicity of notation, where 

r j(0 = r ;(x(0)- 


The derivatives in the global and local descriptions are related by the chain rule 

drj(x) _ / 2 ' \ dTj{ Q 

dx \hj J dt; 


(2.4) 


For simplicity, when the index is not essential, £)■ is denoted by E. Let v and w be two functions on E; 
their inner product is given by 


(v,w) E = I u(x)w(x) dx. 

Je 

The L 2 norm of v on E is (a different norm will be discussed later) 

IMI e = (v(x)) 2 dx 

The inner product in the global description relates to that in the local description by, 

( v , w) E = j u(x)w(x) dx = Jv(0w(0 d% = (^jj ( v , w),. (2.5) 

For any nonnegative integer m , let P m be the space of polynomials of degree m or less. The 
convention P m = (0) for m < 0 is also used. 

Let the Legendre polynomial L m be defined as the unique polynomial of degree m that is orthogonal 
to P m - 1 and L m ( 1) = 1. Then, it is well-known (Flildebrandt 1987, Hesthaven and Warburton 2008) that 
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\\L-m II 2 — (^m< Lm) — ^ 


( 2 . 6 ) 


2.2. Discontinuous Galerkin formulations. Since there are several subtle differences among the 
various standard and new formulations, for ease of comparison, first, the weak and strong forms of the 
DG method are re-derived. 

Let k be a nonnegative integer. At time t, let the solution u(x, t) be approximated on each cell £)■ by a 
polynomial of degree k in x denoted by u ; (x, t). The collection of all {u ; (. , t)} as j varies forms a 
function denoted by u h , which is generally discontinuous across cell interfaces. In the reference 
description for E = £), with t understood, Uj can be expressed in modal form as 

k 

u j {0 = Y J u j.i L ^)- (2-7) 

i = o 

Assume that the data Uj t (t n ) = Uj t are known for all cells. We wish to calculate ^ u ; - ; at t = t n . 

Note first that since u h is generally discontinuous across cell interfaces, f(u h ) defined by f(u h )(x) = 
f(u h (x )) also shares this property. For the data among cells to interact and, simultaneously, the method 
to be conservative, we need to define a unique flux for each interface that is common for the two adjacent 
cells. Also note that, on each cell, f(u h ) can be a non-polynomial function and needs to be approximated 
as will be discussed later. 

Focusing on the cell E = £)■, let 0 be a test function, i.e., an element of P k (independent of t). The 
conservation law (2.1) leads to the following requirement 

— K, 0) £ + ((f(u h )) x , 0)^ = 0. (2.8) 

The above relation on E, however, involves no interaction with the data of the neighboring cells. To 
account for data interaction, first, by integration by parts: 

d 

— {u h ,(l))E + ( f(.u h )(p) dE - = 0. 

For the boundary term ( f(Mh)4>)dE , at each interface j + 1/2, from the left and right u values, namely, 

u ~ = u j+i/2 = u ji x j+ 1 / 2 ) and u + = Uj +1/2 = u j+1 (x j+1/2 ), (2.9) 

a flux value 

f l+i /2 =f I (u~,u + ) (2.10) 

common for the two adjacent cells Ej and Ej +1 is employed (the superscript / represents ‘interface’). This 
common flux, also called numerical flux, is typically an upwind-biased flux for advection-type problems. 
In the case of an advection with constant speed a, 
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( 2 . 11 ) 


1 1 

f I = -a(u + u + ) — ~p\a\(u + — u ) 

where 0 < /? < 1 (with /? = 0 recovering a centered scheme and /? = 1 recovering an upwind scheme). 
For a diffusion-type problem, the common flux is typically a centered quantity. As for 0, the value from 
E (no upwinding) is used. For the DG method, the solution is required to satisfy 

— ( u h , 0) £ + - ( f(u h ), C/) X ) E = 0. (2. 12) 

The above is often referred to as the weak form. Denoting the two boundaries of E by the indices L and R, 

— ( u h , 0) £ + fi - fl 0 L - C f{u h ), (p x ) E = 0. (2.13) 

Let the two boundary fluxes using the data inside the cell E = Ei be 

ft = f(Uj{Xj- 1/2 )) and f R = f(u j (x j+1/2 )). (2.14) 

Again at the two boundaries of E, set 

[/]l = fl ~ fi and [f] R =fk~fR- (2-15) 

That is, [/] represents the jump or difference between the common interface flux and the flux value just 
inside the cell (and [f] L is not f E — f E ). By applying integration by parts again to (2.12), 

— (u h ,0) £ + ((/(u h )) xJ 0)^ + ([/]0)du = 0. (2.16) 

Comparing the above with (2.8), the difference is the boundary term ([f](p)dE- which accounts for data 
interaction. Equivalently, 

— (u h ,0) £ + ((f(Uhj) x .4>) E + [ f]n4>R - [/]l0l = 0. (2.17) 

Note that the weak and strong form s (2.13) and (2.17) are equivalent if the inner products are carried 
out exactly. In practice, however, the inner products are approximated by quadrature rules, and the two 
forms generally produce different results. 

3. New Strong Forms for the DG Method 

The next goal is to rewrite the strong form (2.17) in a manner involving no test function. The 
formulation below deals with f x directly and thus differs from the flux reconstruction (FR) formulation, 
which first reconstructs the flux /. As will be shown, the two formulations yield essentially the same 
schemes with FR typically formulated with the nodal form. The approximate Dirac delta functions below 
also simplify the energy-stability proof in Section 5. 

3.1. Approximate Dirac delta functions. These functions are motivated by the boundary terms in 
(2.17). From here on, unless otherwise stated, inner products take place in the reference frame, i.e., on 
/ = [—1,1]. For a fixed a on /, let the approximate Dirac delta function at a to degree k be denoted by 
8 a k and defined as the linear functional on P k that satisfies, for all 0 in P k , 
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S a ,k (0) = 0(a)- 


(3-1) 


By a theorem in functional analysis, there exists a member y a ,k ' n Pk such that for all 0 in P k , 

(Ka,fe-0) = 5 a , fe (0) = 0(a)- 

In fact, Ya,k can t> e calculated. First, since it is of degree k, it can be expressed as 

k 

y«, fe (0 = ^\Li(0- 

;=o 


(3.2) 


(3.3) 


For any m < k, to obtain b m , apply the dot product with L m to the above and, using the fact that L m is 
orthogonal to all L;, i ^ m, 


(ja.k'Lm) ~~ 


By (3.2), ( Ya,k’ Lm) = L m (a). This fact together with (2.6) imply 


^m(a) 


2 

2m + 1 ' 


That is, b m = — - — L m (a). Consequently, 


L t (f). (3.4) 

i= 0 

See e.g., (Morse and Feshbach 1953, p. 719) where an expression for the Dirac delta function in terms of 
general eigenvectors can be found. 

For simplicity, S a k is considered identical to y a , 



a ^a,k Y a,k Y a- 


Note that fora = +1, |L;(a)| = 1; therefore, L;(a) 2 = 1, and 


A. 

i,fc(i) = ^-i,fe( _ i) = 'y' 

i = 0 


2i + 1 
2 


As a result, 


SiAV = s-i, fc(-i) 


(k + l) 2 
2 


(3.5) 
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The above is the maximum value of <5 ±1 on / and is related to the fact that in conjunction with an explicit 
time-stepping method, the CFL condition for the DG method is inversely proportional to ( k + l) 2 . In 
addition, the i-th mode of 5 ±1 satisfies, by (2.6), 

2i + 1 

— 2~ ii(± !) k 

Therefore, after some algebra, 

(3-6) 

The graphs of 8 0 k , k = 0, 2,4,6, and 8 are shown in Figure 3.1(a). Note that since L 2k+1 ( < 0) = 0, 
^o, 2 fc+i = ^ 0 , 2 k- F° r comparison, the graphs of 5 ljfc , k = 4, ...,8 are shown in Fig. 3.1(b). Note the 
difference in scales. 

3.2. Strong form SI. We can now derive the new strong forms. Again unless otherwise stated, we use 
the local description, and the indices L and R represent —1 and 1, respectively. By (2.4) and (2.5), the 
strong form (2.17) implies 



yj (%,0) t + ((/(u/i)),, 0) + [f] R (p R ~ [f] L 4>L = 0. 


(3.7) 


Using (3.4), set 


k k 

Z 2i + 1 v -1 2i + 1 

— 2 — (-1 yL ^° and = 


(3.8a, b) 


m= 0 


m= 0 


Approx. Dirac delta function at x = 0 

y 


Approx. Dirac delta function at x = 1 

y 



k = 4, 5,. ..,8. 
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Then S L and S R are of degree k and, for all <p in P k , 


(5 l ,0) = 0 l and (S R ,(p) = <p R . (3.9) 

Thus, (3.7) and the above imply 

(jj {u h ,<t>\ + ((f(u h )) f , 0) + [fl R (S R ,<P) - [f] L (S L ,<P) = 0. (3.10) 

For any integer m > 0, denote by T m the projection onto P m , i.e., with any integrable function v, 

m . m 

Z (v,Lj ) v 1 2i + 1 

L i =}_ l — — ("> L d L i- (3.H) 

i = 0 ‘ i = 0 

Due to the fact that 0 is in P k , (3.10) is equivalent to 

(j) (u h ,4>) t + (^ k [(/K)) f ],0) + [/] r (5 r ,0) - [f\ L (S L ,<p) = 0 . (3.12) 

The projection T k above serves the purpose of cancelling out the test function <p. Indeed, since u h , 
^fc[(/( u h))f]’ ^r- an d f>L ar e all in P k and, since <p spans P k , the above implies the following strong 
form referred to as SI, 

SI (y) (u h )t + ^ k [(/(%))f] + [HrSr ~ [f] L S L = 0- (3-13) 

Two remarks concerning the form SI are in order. 

1. Generally, after integrating by parts twice, we get back the original equation. Here, due to the 
different flux values employed for the boundary terms in each of the two integrations by parts, the strong 
form SI looks similar to but not identical with the original conservation laws in the local description. The 
differences are: (a) the derivative with no interaction is provided by Pk[(f( u h))$] and not just ( f(u h ))^ 
and (b) the correction for the derivative, namely, \f\ R fi R — \f\ L fi L , provides data interaction among cells 
and, at the same time, deals with the flux jumps at the interfaces. 

2. For numerical algorithms, the interaction of the data among cells involves the functions S L and S Rl 
which are derived exactly and are given by (3.8). Concerning P k [(f(u h ))^], which represents the 
derivative with no interaction, u h is in P kl but f(u h ) needs not be a polynomial in general. For nodal- 
type schemes, P k [(f(u h ))^] can be approximated conveniently and effectively by using the chain rule at 
the nodal points ft, i = 0, ...,/e: (/(u h )) (ft) = a(u h (ft)) (u h ) f (ft). 

3. The strong form (3.13) leads to a family of schemes by altering the criteria of approximation for the 
Dirac delta function. What is crucial is that the resulting scheme must be stable and accurate. 

3.3. Strong form S2. We next present an alternative to the strong form SI. Concerning again the 
derivative with no interaction, we can first approximate /(n h ) by its projection P k (f(Uh)) and then 
calculate the derivative: the result is referred to as S2, 
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S2 


yl ( u n)t + [^fe(/( u h)]f + IHr^r ~ WlSl ~ 0- 


(3.14) 


For this strong form, the two boundary fluxes interior to cell E = £)• are defined as, instead of (2.14), 

/l + = [?k (/(«/))] (*/- 1 / 2 ) and = £P k (/(uy))](xy + i/ 2 ) (3.15) 

and the jumps [/] at the two interfaces are defined using the above interior values. 

The following remarks are in order. 

1. For this strong form, similar to the case of SI, concerning the derivative with no interaction 
Pfc(/( u /i)]f> with nodal-type schemes, we can first approximate f(.u h ) by its values at the nodal points 
and then evaluate the derivative of the resulting polynomial. Such an approximation is said to be via 
‘Lagrange polynomial’ whereas that in the second remark after (3.13) is via ‘chain rule’. These two ways 
of approximating the derivative with no interaction were mentioned by the author as alternative 
approaches (Huynh 2007). Later, Wang and Gao (2009) found that for nonlinear problems such as the 
Euler equations, the solution via the chain rule was more accurate. 

2. With the current strong forms SI and S2, the conservation law for nodal-type methods is solved via 
a combination of interpolation and projection: the derivative with no interaction is approximated via 
interpolation, and the two boundary corrections are carried out exactly to degree k via projection since the 
approximate Dirac delta functions <5 ±1 are exact up to degree k. The evaluation of <5+i(0 at the nodes 
incurs no additional (interpolation) errors. Also note that the approximation to the derivative with no 
interaction (/(n ft ))^ can be carried out via projection as well but, typically, projection for such a function 

is costlier than interpolation; moreover, projection is typically approximated by a quadrature which, 
similar to interpolation, incurs errors. 

3. For the case of advection with constant speed, the strong form s SI and S2 become identical and is 
called strong form S. The reason is that since /(u.) = auj is of degree k, both approaches for the 
derivative with no interaction result in a(u ; )^, which is of degree k — 1. In the general case of nonlinear 
conservation laws, however, the two forms are different: the polynomial \'E’ k (f (u h )\^ for S2 is of degree 
k — 1 as opposed to Pk[(f( u h))%] for SI, which is of degree k. The lower degree approximation of S2 
appears to be consistent with the numerical findings of Wang and Gao (2009) mentioned above. 

4. The key difference between the current strong forms and a standard DG strong form such as (2.17) 
is that by employing the approximate Dirac delta functions, the current ones involve no test functions. As 
will be shown, they also unify existing methods by recovering the spectral difference schemes via an 
appropriate modification to the approximate Dirac delta function; in addition, these strong forms result in 
a family of energy-stable methods that have the property of super accuracy (i.e., more accurate than 
expected via Fourier analysis as discussed by Huynh (2007) and Vincent et al. (2011)). 

3.4. Conservation property. Next, we show that the strong form SI is conservative, i.e., 

^ J u h (x) dx = fl - fl (3.16) 

Note first that for any integrable function w on /, the projection preserves the cell average quantity, 
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By (3.23) and (2.15), Eq. (3.20) follows. This completes the proof. 

The conservative property of the scheme S2 is similar with (3.15) replacing (2.14) and, therefore, 
omitted. 

A modification to assure conservation property when using the derivative with no interaction via the 
chain rule on 2D unstructured meshes was provided by Gao and Wang (2013). 


4. DG Scheme Via Flux Reconstruction (FR) Approach 

We now present a formulation which deals with / rather than f x as in the previous section. It will be 
shown that for the DG method, the slope f x for the conservation law is given by (Fj) x where, on each cell 
Ej , Fj is a polynomial of degree k + 1 defined by (a) Fj takes on the common interface flux values at the 
two cell boundaries and (b) T k _ 1 {Fj) = F k _ 1 (f(u h )'). The key idea here for constructing Fj is to 
integrate the strong form S 1 . 


4.1. Approximating the flux /(u h ) with no interaction. In the form SI, the derivative with no 
interaction is approximated by !P fe [(/(u h ))^]. This approximation is equivalent to a certain approximation 
for the flux / itself. In what manner is / approximated? 

To answer the above question, first, by the fundamental theorem of Calculus, 

(/(tfiJX??) = fh +j (f(u h )) ( df. ( 4 . 1 ) 

Next, with subscript IPD short for ‘integrating the projection of the derivative’, set 

/ ipd 07 )=/ l + + [ P k [(f(u h ))t] ( 4 . 2 ) 

J -i 
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Then /j PD is of degree k + 1 and, the second equation below is due to (3.19), 

/ipd( - 1) = fi an d /ipd(1) = /r • (4-3) 

That is, / IPD inteipolates f(u h ) at the two cell boundaries. Two conditions for / IPD are known, k 
conditions remain to describe how it approximates f(u h ). 

To obtain these remaining conditions, let <p be in P k . Integration by parts yields 

(/ipd. ) = /r 0r — fi 0l ~ ((/lPD)f. 0)- (4-4) 

By (4.2), (/ IPD ) C = TkKfiUh))^]. Thus, the above implies 

(/ipd- 0f) = /r" 0r - fh (Pl ~ (^[(/OfiJ)?], 0)- (4.5) 

Next, in (4.5), since 0 is in P k , we can replace ^[(/(Wi))^ by (/(u h ))^, and since 0^ is in P k - 1} we 
can replace / IPD by ^ fc -i(/i PD ), 

(^ViC/ipd), 0f ) = /r" 0r - A + 0l - ((/Oh))f . 0)- (4-6) 

By applying integration by parts to the inner product of the right hand side, 

C^-iG/pd). 0f) = ((/(%)), 0f). (4.7) 

As 0 spans P k , 0^ spans P k -\, therefore, the above implies 

^ fe-l(/lPD) = P fc-l(/( u h))- (4-8) 

Thus, in the strong form SI, T^l (/(iq,)) | approximates ( f(u h ))^ in the sense of projection by its 

definition. Equivalently, its integral / IPD defined by (4.2), which is a polynomial of degree k + 1, 
approximates f(u h ) in the sense of a combination of interpolation at the two boundaries as in (4.3), and 
projection as above. 


4.2. Correction functions. We now deal with the boundary terms of the form SI or (3.13). With S R by 
(3.8b), let the correction function for the right boundary be denoted by g Ri DG or simply g R and defined by 



S R (g)dg. 


(4.9) 


(In Section 6, g R denotes the correction function for a more general scheme, not just DG.) Then, g R is of 
degree k + 1 and 


9r ~ Sr- 


(4.10) 


Moreover, 


g R (~ 1) = 0 and g R { 1) = 1, 


(4.1 la, b) 
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where the second equation results from (4.9) with ^ = 1 together with (3.2) where 0 = 1. Note that since 
8 r approximates the Dirac delta function at ^ = 1, g R approximates the Heaviside or unit step (step up) 
function at the same location. 

Next, it is shown that 


= 0. (4.12) 

Indeed, by integration by parts, for any polynomial v of degree k, 

1 1 

(g R ,V) = KD - [ gktfMOdf = u(l) - [ = 0 . (4.13) 

4-i J-i 

As v spans P k , v' spans P k _ therefore, (4.12) holds. 

Let the correction function for the left boundary be defined by (note the change in the integration 
limits) 


gdO = f 5 l 0?)*?- (4.14) 

(Again a more precise notation is ,g L DG ; however, the subscript DG is understood in this section.) Then, 
g L is of degree k + 1 and 


gL d L . 


(4.15) 


Moreover, 


0i(-l) = 1 and g L ( 1) = 0, (4.16a,b) 

where the first equation results from (4.14) with <f = — 1 together with (3.2) where 0 = 1. Note that since 
8 l approximates the Dirac delta function at £ = — 1, g L approximates the step down function which takes 
on the value 1 at ^ = —1 and the value 0 for —1 < ^ < 1. 

In addition, 


^ fc-i( Sf.) — o. 

The proof is similar to that of (4. 12) and is omitted. 


(4.17) 


4.3. Flux reconstruction. The reconstruction of the flux on E = Ei can now be carried out. At the two 
interfaces, recall that / L + = /(uy(-l)), f R = f(uj( 1)), [f] L = f[ - / L + , and [f] R = f R - f R . 

From Subsection 4. 1 above, / IPD is of degree k + 1, involves no interaction among cells, takes on the 
values fi and f R at the left and right interfaces, respectively, and T’ fe _ 1 (/ IPD ) = P k -i(f(u h )). Next, the 
function 


/ipd + \f] L gL 
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recovers the common flux value at the left interface, namely fj_ 1 / 2 , due to the condition g L (—l) = 1. It’s 
value at the right interface is unchanged from that of / 1PD , namely the value f R ; this is due to the 
condition g L ( 1) = 0. The function 


Fj = fiPD + if] L 9l + [f] R 9r (4. 1 8) 

takes on the common interface flux values at the two boundaries, 

Fj{-l)=fl and Fj(l)=fl (4.19) 

In addition, on Ej, by (4.8), (4.12), and (4.17), 

? k -i{Fj) = ^(/(u,)). (4.20) 

Finally, by (4.2), (4.10), and (4.15), Eq. (4.18) implies 

Fj' = r k [(f(Mh))f] - WlSl + [fhS R . (4.21) 

Note that the above right hand side is identical to that of the strong form SI of (3.13). This completes the 
derivation of DG schemes via FR approach. 

A few remarks are in order. 

1. Due to (4.19), the collection of all {Fj} as j varies forms a function which is continuous across cell 
interfaces, and Fj' yields the DG approximation for f x for the conservation law. 

2. The flux reconstruction form (4.18) above leads to a family of schemes by using different criteria to 
define the correction function g L and g R . What is crucial is that the criteria must result in stable and 
accurate methods. 

3. In the weak form (2.13) for the DG scheme, f(u h ) can be replaced by since <p x is of 

degree k — 1, 


a 

— (%, 0)u + fl cp R - fl 0 L - CPfc-1 (/(%)), 0 X ) £ = 0. (4.22) 

Thus, due to (4.19) and (4.20), the definition for Fj by (4.18) is, loosely put, consistent with the weak 
form above. 


5. Energy Stable Flux Reconstruction Schemes (ESFR) 

Energy stability for a family of schemes can now be shown. Consider the case of advection with 
constant speed a > 0 (the wind blows from left to right). We also assume that the initial data is periodic 
or of compact support so that boundary conditions are trivial. 


5.1. Energy stability of the DG scheme. First, it is well-known that the DG scheme with the upwind 
biased flux (2.11) is energy stable. The simple proof, which can be found in (Jameson 2011), plays a 
critical role below; therefore, it is included. 

The strong form (2.17) on E = Ej with 0 = iq implies 
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fj-i/ 2)u/-i/ 2 } = 0. (5.1) 


(^t Uj ’ Uj ) E + ^ ( Uj '^x’ U j) E + (fj+l/2 ~ fj+l/2) U j+l/2 ~ Uj-t/2 ~ 

Note that (a(Uj)f, Uj) = ^(ii“ +1 y 2 ) 2 — ( u /-i/ 2 ) 2 ]- With /y+1/2 = au j+i/2 f° r eac h interface, the 
expression in the curly bracket above can be written as 

a ( u j + 1/2 — 2 U j + l/2) U j+l/2 ~ a ( u j- 1/2 — 2 u j-i/ 2 ^ u j-i/2- (5-2) 

At each interface 7 — 1/2, collecting the contributions from the element Ej_ 1 and then Ej, there is a total 
contribution of 


a (Uj_i/ 2 2 U j-l/2) U j-l/2 a (Uj_ 1/2 2 u f-l/2) U f-l/2- (5-3) 

1 1 

Omitting the subscript j — 1/2, we obtain a[(u 7 — -u~)u~ — ( u 1 — -u + )u + ]. After some algebra, and 
using the upwind-biased flux (2. 1 1), the result is 

^/3\a\(u + -u~) 2 . (5.4) 


Recall that 0 < /? < 1. Thus, on the whole physical domain, with our periodic or of compact support 
assumption for the data, (5.1) and the above imply 


1 d 

2 at 


1 


ut 


dx = — 


|a|(u/_ 1/2 


U j — l/2 



(5.5) 


That is, energy of the numerical solution can only decrease in time. This completes the proof that the DG 
scheme in semi-discrete form using the upwind-biased common flux is energy stable for the linear 
advection equation. 


What we need from the above proof is the fact that the quantity in the square bracket of (5.1) implies 
(5.4) at each interface after the square bracket quantities from the two adjacent cells are collected. 


5.2. Energy stability of a family of FR scheme. For the current case of advection, as discussed in the 
third remark after (3.15), the two strong forms SI and S2 become identical and are referred to as strong 
form S. On Ej, set 


fj = f(uj ) = auj. (5.6) 

To obtain a family of schemes, we replace 8 in the strong form S of (3.13) by 8 + e, where e provides the 
variation: 



( u j)t + [/]l(<5l + e i) + [/]r(<5r + e R ) - 0. 


(5.7) 


In other words, we approximate the Dirac delta function not by 8 but by 8 + e, which results in a family 
of schemes. What is crucial is that e must be chosen in a manner that the resulting scheme is stable and 
accurate. In this section, we assume e involves only the highest mode: 
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and e R = e k L k . 


(5.8a,b) 


e L — ( ^) ke k^k 

The coefficient (— l) fe above is for being consistent with (3.8a). We wish to show that if 


Eb- > 


2/c + l 


(5.9) 


then the resulting scheme is energy-stable under a certain energy > norm to be defined later. 

2/c + l 

Compared with the coefficient of the highest mode of S R which, by (3.8b), equals — — , the above 
condition is equivalent to the fact that the highest mode of S R + e R has a strictly positive coefficient, i.e., 


k- 1 


Z 2i + 1 

— 2 — L-i = a k Lk + S R , k - i ■ 
l = o 


(5.10) 


where 


2/c + l 

a k = 2 + £k > °- 


(5.11) 


For the left boundary, 


$l + e L = (-1 Ycc k L k + . 


(5.12) 


The following proof is similar to but simpler than that of Vincent et al. (2011) due to the use of the 
approximate Dirac delta functions. 

Using technique employed in (5.1), multiplying (5.7) by Uj and integrating on /, we obtain 


Set 


hj d 
4 dt 



~L u 4fj) ( 


df + lf] L u, iL - If] 


R U j,R 


+ 




Uje R d^ . 


(5.13) 


R = 




Uj E r . 


(5.14) 


Eq. (5.13) implies 

~4 u f d % ~ R = ~\^j u j(fj){: di > ~ IHlUj-.l + [f] R Uj, R J. (5.15) 

The square bracket quantity above equals the square bracket quantity in (5.1), which will result in (5.4) at 
the interfaces. 
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We now derive conditions on e k so that the above left hand side results in the time rate of change of a 
certain energy norm. Our goal is to relate R of (5.14), which does not involve with 7 ( Uy fe (t)) where 
uj k is, up to a constant factor, the energy of the k- th mode of Uj. 

To this end, first, on /, 

^ k k 

Hl 2 =f u f d t = = ^2TTT U7?; ' (5-16) 

-1 i=0 i = 0 


Next, since e L and e R involve only the highest mode as expressed in (5.8), and since the k- th mode of 
Uj is Uj k L k , Eq. (5.14) implies 

R= 6 k u jik \\L k \\ 2 {(-l) k [f] L -[f] R l (5.17) 

Following (Jameson 2010) and (Vincent et al. 2011), differentiate (5.7) k times (in space), 


hj d d k Uj 


) /c + 1 


+ ■ 


fj 


2 dt df k 


fc + 1 


- [/] 


d k (S L + e L ) 
df k 


+ ifl 


d k (S R + e R ) 
df k 


= 0. 


(5.18) 


Since fj is of degree k. 


Multiplying (5. 18) by 

hj d fd k Uj \ 2 

7 at \ 7 J* ) 

Integrating on /, we obtain 


where 


d k+1 fj 

df k+1 


Ifi 


d k uj d k (S L + e L ) 


df k 


df k 


+ [/L 


d k uj d k (S R +e R ) 


df k 


df k 


= 0. 


hj d r 1 (d k Uj\ 2 

7 It J_! j 


= R 


(5.19) 


(5.20) 


(5-21) 


R = 



d k (S L +e L ) 
df k 


df 



d k (S R +e R ) 
d% k 


df. 


(5.22) 


For the left hand side of (5.21), denote 

d k L k 


Then, m k is a constant. In fact, m k = k\ a k = (2k — l)!/(2 fc 1 (/c — 1)!) where a k is the leading 
coefficient of L k ; however, this expression is not essential. Concerning the left hand side of (5.21), 
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(5.23) 


d k Uj d k L k 

g^k = U j,k g^k = U j,k m k- 

Consequently, recalling that Uj k = u ; - fe (t), the left hand side of (5.21) equals y y (m k uj k ) where we 

f 1 

use the fact that J 1 = 2. That is, (5.21) can be written as 


hj d r 1 /d fe u ; \ 2 


df 


hj d 
2 dt 


ulk) = R. 


(5.24) 


We now relate R defined by (5.22) with R of (5.17). The key fact employed is that R involves only the 
highest mode. First, by (3.8a) for S L and (5.8a) for e L , 


d k (S L + £ L ) 

d% k 


, (2k + 1 \ 

= (_1) 


m k> 


and by (3.8b) for S R and (5.8b) for e R , 

d k (S R + e R ) (2k + 1 


(2k + 1 \ 

d? ={-^- + ek ) mk - 
Due to the above two expressions and (5.23), Eq. (5.22) implies 

R= (2/c + 1 + 2e k ) m 2 k u j k { (-1 ) k [f] L - [f] R } 
Using the above and the fact that ||Lfc|| 2 = 2/(2 k + 1), Eq. (5.17) implies 

2e k 


R = 


(2k + 1)(2 k + 1 + 2 e k )ml 


R. 


Thus, by (5.24), 


R = 


h, 


d(Uj, k ) 


(2k + 1) (2/c + 1 + 2e k ) ’ dt 


(5.25) 


(5.26) 


(5.27) 


What is crucial is that R of (5.14), which involves no time derivative, can be expressed as above, which 
involves the time derivative and the highest component Uj k . 

By the above, the left hand side of (5.15) yields 


hj d 
4 dt 



h, d 

n — _L 

4 dt 



46fc 

(2k + 1) (2/c + 1 + 2e k ) 



(5.28) 


Using (5.16), the k - th mode of the above right hand side equals 


hj d 
4 dt 



2 ffc 

(2k + 1 + 2e k ) 


2k + 




For the above to result in a norm, it is required that 


(5.29) 
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(5.30) 


1 - > 0 

(2k + 1 + 2e k ) 

or, (2 k + 1 + 2e k — 2e k )/(2 k + 1 + 2e k ) > 0 or (2k + 1) / (2 k + 1 + 2e k ) > 0 or (2k + 1 + 2e k ) > 
0, or 


e k > 


2k + 1 
2 


(5.31) 


Thus, instead of the standard norm (5.16), we employ a modified energy norm on E = £)■ defined by 


u f = 

7 ll M 




26b 


(2k + 1 + 2e k )J 2k + 1 


u j,k ~f 


k- 1 
i = o 


— 7 U / i 
+ 1 ; ' 1 


(5.32) 


Or, equivalently, by (5.23), 


U; = = 
7II m 


J uf cL% — 


26b 


(2k + 1 + 2e k )(2k + 1 )m 


i: 


d k u } 

~d^ 




(5.33) 


Thus, using the above norm, the family of schemes via the approximate Dirac delta function of <5 + 6, 
where the highest mode of S R + e R = a k L k + S R k _ 1 has a strictly positive coefficient a k , is energy 
stable. This completes the energy-stability proof. 

The following three remarks are in order. 

1. The energy-stability proof of Vincent et al. (2011) was extended to nonlinear conservation laws in 
Jameson, Vincent, and Castonguay (2012). There, it was also shown that similar to DG, the FR methods 
may encounter aliasing driven instability if the flux function is non-linear. They demonstrated that for 
nodal methods, the location of the solution points plays an important role in determining the extent of 
aliasing driven instabilities, and the strong quadrature points such as Gauss points are optimal in 
minimizing such instabilities. 

2. The energy-stability proof here also extends to nonlinear conservation laws provided that (a) the 
strong form S2 is employed and (b) the interface flux is an E-flux (Osher 1984, Jiang and Shu 1994); this 
extension is similar to that by Jameson, Vincent, and Castonguay (2012) mentioned above. 

3. It remains an open problem to prove or disprove that for the strong form SI, in the case of nonlinear 
conservation laws, the family of schemes discussed here is stable. The above proof does not extend for 
the following reason. A key ingredient of the proof is condition (5.19) or d k+1 fj/d^ k+1 = 0. For the 
strong form SI of (3.13), the quantity d k+1 fj/d^ k+1 is replaced by d k (T k [(f (u h )) ^]) / d^ k , which may 
not vanish. Also note that for nonlinear conservation laws and the strong form S 1 , whereas the stability of 
the family of schemes is not known, the stability for DG using the weak form (and thus the strong form 
SI) is well established as discussed in (Jiang and Shu 1994). 


6. Flux Reconstruction/Correction Procedure via Reconstruction 

(FR/CPR) Schemes 

The FR methods, which are obtained by different criteria of approximation for the (Dirac) delta 
function in the strong form SI, or different definitions for the correction functions in constructing f) via 
(4.18), are discussed in this section. Their accuracy and stability properties by Fourier analyses are also 
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mentioned. The current presentation of these methods provides additional relations not available in the 
literature. 

Focusing on the right boundary, recall that in the case of the DG scheme, the correction function 
9r,dg is °f degree k + 1 and defined by the k + 2 conditions of g R DG (— 1) = O' S'r.dgCI) = 1' and 
Pk-iidR.Dc) = 0 of (4.12). Its derivative is the approximate delta function S R = S R k given in (3.8b). 

To obtain a family of schemes, the correction function g R of degree k + 1 is required to satisfy 

g R (-l) = 0, g R ( 1) = 1, (6.1a,b) 

together with k additional conditions. As discussed before (4.18), when reconstructing the flux, the 
condition g R ( 1) = 1 serves the purpose of recovering the common flux value fj+1/2 at the right interface 
whereas condition g R (— 1) = 0 leaves the value at the left interface unchanged. The k additional 
conditions correspond to the fact that g R approximates the zero function on [—1, 1) (consistent with the 
remark after (4.11) that g R approximates the unit step function) in some sense. These approximation 
criteria must be chosen in a manner that the corresponding scheme is stable and accurate. 


6.1. Radau and Lobatto polynomials and their derivatives. The following polynomials play an 
important role in the definitions of various correction functions and their derivatives. 

Radau polynomials. Let the left Radau polynomial of degree k + 1 be defined by 

R L , k+1 =^L k+ i + L k ). (6.2) 

Whereas L k+1 vanishes at the k + 1 Gauss points, R L k+1 vanishes at the k + 1 left Radau points. At the 
boundaries, 


R 


L,k+1 


(-1) = 0 and R L , fe+1 (l) = l. 


(6.3a,b) 


Since both L k and L k+1 are orthogonal to P fc _ 1 , Ri , k + 1 also possesses this property. Equivalently, 

^ fe -iOW) = 0. (6.4) 

Equations (6.3) and (6.4) above imply that R L k+1 is the correction function g R DG . Thus, by (4.10) and 
(3-8), 

k 


R 


L, k + 1 


' - S Rjk - ^ 


2i + 1 


L, 


m = 0 


In a similar manner, the right Radau polynomial of degree k + 1 is defined by 

(_ 1 ) fe+ 1 


R 


R, fc + 1 — 


'(ife+i — Lk)- 


(6.5) 


( 6 . 6 ) 


Then R r k+1 vanishes at the k + 1 right Radau points. In addition, 


R 


R, fc + 1 


(-1) = 1, R R , k+1 ( 1) = 0, 


(6.7a,b) 


and 
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^-i(^, fe+ i) = 0. (6.8) 

Properties (6.7) and (6.8) above imply that R Rk + 1 is the correction function g L DG . By (4.15) and (3.8), 


R 


R,k + 1 


'=-^=-X 


2i + 1 


(-1 )% 


(6.9) 


m= 0 


The plots of the first few left and right Radau polynomials and their derivatives are shown in 
Figures 6.1 and 6.2. 
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Figure 6.1. — (a) Graphs of the left Radau polynomials R L k , k< 4 (the dots are the 4 left Radau points) and 
(b) their derivatives. The k - 1 zeros of R Lk together with x = 1 form the k right Radau points. 
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Derivatives of Right Radau Polynomials 
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Figure 6.2. — Graphs of (a) the right Radau polynomials R r k , k < 4 (the dots are the 4 right Radau points) and 
(b) their derivatives. The k - 1 zeros of R Rk together with x - -1 form the k left Radau points. 


Lobatto polynomials. For any integer k > 2, let the Lobatto polynomial of degree k be defined by 
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( 6 . 10 ) 


Lc>fc — L k — L k _ 2 - 


Then 


Lo fc (-l) = Lo k (l) = 0. (6.11) 

and, by definition (6.4), Lo fe is orthogonal to P fc _ 3 . The zeros of the Lobatto polynomial of degree k are 
the k Lobatto points; they include the two boundaries +1 . The derivative of Lo fe+1 , which has not been 
employed by FR authors and play an important role here, is given by 

LOfc+i/ = (2k + V)L k . (6.12) 

The plots of the Lobatto polynomials Lo fc , k = 1, ... , 4, and their derivatives are shown in Figure 6.3. 

6.2. FR schemes. Due to symmetry, we focus only on g R (g L involves an extra factor (— 1)*). For the 
DG scheme, as discussed after (6.14), 

9r. DG = ^L,k+l- (6.13) 

Note that the correction function for the right boundary is the left Radau polynomial. 

The correction functions for the ESFR methods of Section 5 can now be derived. First, for simplicity, 
denote 


9 ~ 9r- 


Lobatto Polynomials 


Derivatives of Lobatto Polynomials 
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Next, g is assumed to satisfy the two conditions (6.1). For the k additional conditions, instead of 
requiring g to be orthogonal to P k -\ as in the case of the DG scheme, we require it to be orthogonal to 
P k _ 2 . It was verified via Fourier analysis (Fluynh 2007, 2009) that if g is orthogonal to P k - 2 the 
resulting scheme is stable (the converse is not true, however). 

Since both R Lik+ i and R Lk are orthogonal to P k ~ 2 , if g orthogonal to P fc _ 2 , it can be written as 

9 = aR Ll fc+i + (1 - a)R L , k (6. 14) 

where a remains to be specified (note that the parameter a here is completely independent from the a for 
the location of the approximate delta function in Section 3). The current author considered only 0 < a < 
1 in (Fluynh 2009) for the obvious reason that such a correction function is a weighted average of the two 
Radau polynomials (correction functions). For a > 1, the correction function becomes steeper than g DG 
and, in conjunction with an explicit time stepping method, the largest possible time step size or CFL limit 
becomes smaller than that for DG, a disadvantage. For a < 0 the scheme is unstable. Concerning only 
energy stability and accuracy, however, a can be larger than 1. 

The above implies, 


r \ a 

g = cc{R Lk+ 1 — Ri :k ) + R k)k = — Lo k+1 + R Rk . (6.15) 

The second equality is a result of (6.10) and (6.2). Note that a = 1 recovers the DG correction function. 
By (6.12), the above implies 


g' = |(2 k + 1 )L k + R[, k = | (2 k + 1 )L k + S Rik . v 


Thus, by (5.10), we obtain the ESFR methods with 


a 


a k = -(2 k + 1), 


(6.16) 


(6.17) 


and the method is energy stable if a > 0. 

Concerning accuracy, it was found by experiments using Fourier analyses (Fluynh 2007) that if g is 
orthogonal to P m , then resulting scheme is accurate to order k + m + 2. Using this criteria, the DG 
scheme is accurate to order 2k + 1, and the schemes of the ESFR family, order 2k. The expected order of 
accuracy is k + 1. Due to this higher than expected accuracy, these methods are said to be super-accurate 
or to possess the property of super-convergence. 

Next, noteworthy members of the one parameter family of correction functions (6.14) are discussed. 

The first choice for g is (6.13), which results in the DG scheme. 

The second choice for g , denoted by g 2 or c/Lump, l<> (for Tumping for Lobatto points’), is defined as 
follows. Since a steeper correction function tends to result in a scheme with a smaller CEL limit, to make 
g less steep, the extra condition is obtained by pushing one of the zeros to the left boundary, i.e., <f = — 1 
is a zero of multiplicity two. After some algebra (Fluynh 2007), 


9 2 ~ STump, Lo 


l: ^ fc + 1 n 

2fc + 1 Rl ’ k+1 + 2fc + 1 Rl ’ k ' 
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The function g 2 has the following remarkable property. Among the k + 1 Lobatto points (and we can 
use them as solution points), g' 2 vanishes at k of them; the exception is the right boundary where the 
right flux jump locates. 

The final choice here for g requires that g vanishes at the k Gauss points (Huynh2007), 

k + 1 n k 

3Ga = 2k + i Rl ’ k+1 + 2k + 1 Rl ' k ' 

Whereas the staggered-grid (Kopriva and Kolias 1996) and SD schemes are mildly unstable, the above 
provides a modification using the k Gauss points as flux points. The resulting scheme is both Fourier and 
energy stable for all k. 

Also note that the steepest slope of g, which takes place at the right boundary, relates to the time step 
size or CFL limit. For example, for DG, g' = R' L k+ i(l) = (k + l) 2 /2; when an explicit Runge-Kutta 
method is employed, it is well known that the CFL limit for the DG scheme of degree k is roughly 
proportional to 1/A: 2 . 

As an example of correction functions, Fig. 6.4(a) shows the plots of the functions g DG , g Ga , and g 2 
for k = 3. Figure 6.4(b) shows the (top half of) spectra of the corresponding schemes. Note that all three 
methods are stable since the spectra lie in the left half of the complex plane. The spectra intersect the real 
axis at x DG= — 19.2, x Ga= — 12.3, and x 2= — 9.6. If the RK4 method is employed for time stepping, then 
the CFL limits are approximately .145, .227, and .289, respectively. 

Table 1 tabulates the orders of accuracy and errors of these semi-discrete schemes for k = 3. Here, the 
coarse mesh error corresponds to w = n/A, and fine mesh error, w = n/8 so that errors are away from 
machine zero, and the calculation of order of accuracy is valid. 

Note that for each k, all errors have negative real parts — a fact consistent with the stability of the three 
schemes. In addition, the dominant part of the error for the DG scheme is real (dissipation type) whereas 
those for the other two schemes are imaginary (dispersion type). 



(a) Correction functions 


Im 



Figure 6.4. — (a) Correction functions for k = 3, and (b) their spectra. 
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TABLE 1.— ORDERS OF ACCURACY AND ERRORS OF SCHEMES FOR k = 3 


Scheme 

Order of 
accuracy 

Coarse mesh error, w = nl 4 

Fine mesh error, w=; r/8 

DG 

7 

-1. xlO -7 -1.x 1 0 — 8 / 

-4. xlCT 10 -2. xlO -11 ; 

t?Ga 

6 

-3.1 xlO -7 +1.3xl0“ 6 / 

-1.2 xlO -9 + 1.1 x 10 -8 / 

&2 

6 

-5.4x10“ 7 +2.3x10“ 6 / 

-2.2xl0“ 9 +1.9xl0“ 8 i 


A brief review of Fourier analysis for these type of methods can be found in the Appendix. 

Concerning open problems, it remains to be shown whether Fourier stability is equivalent to energy 
stability. 

Also note that there exists correction functions that are not orthogonal to any P m , but the resulting 
scheme is Fourier stable. As an example, the correction function g R that lumps all corrections to the right 
boundary when the solution points are the Chebyshev-Lobatto points (Fluynh 2007) possesses this 
property. It remains an open problem to identify all correction functions that yields stable schemes. 

7. Conclusions and Discussion 

In conclusion, a formulation for the discontinuous Galerkin (DG) method that leads to solutions using 
the differential form of the equation was presented. The formulation includes (a) a derivative calculation 
that involves only data within each cell with no data interaction among cells, and (b) for each cell, 
corrections to this derivative that deal with the jumps in fluxes at the cell boundaries and allow data 
across cells to interact. The derivative with no interaction is obtained by projection or interpolation, and 
the corrections are derived using the approximate (Dirac) delta functions. The formulation results in a 
family of schemes: different approximate delta functions give rise to different methods. It was shown that 
the current formulation is equivalent to the flux reconstruction (FR) formulation. Due to the use of 
approximate delta functions, an energy stability proof simpler than that of Vincent, Castonguay, and 
Jameson (2011) for a family of schemes was derived. Accuracy and stability of resulting schemes were 
examined via Fourier analyses. Some open problems were discussed. The current and FR formulations 
appear to be promising toward the goal of faster and more accurate solutions of conservation laws for 
practical applications. 
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Appendix — A Brief Review of Fourier Stability and Accuracy Analysis 


Accuracy and stability of a numerical method are critical and can be derived via Fourier analyses. On 
(— 00 , oo), consider the advection equation 


u t + u x = 0- (A.l) 

The initial condition at t = 0 is periodic: u init = e lwx where w is a real number between — n and n called 
a wave number and, in this section, / = V— 1. Low frequency data corresponds to w of small magnitude, 
high frequency, to w near +7T. The exact solution is u exact (x, t ) = At x = 0 and t = 0, 

(“exact) t (0,0) = —Iw. (A.2) 

The cells are £) = [/' — 1/2, j + 1/2]. At time t = 0, the data u init on £) is e /w O + ?/2) where |^| < 1. 
Thus, the exact data on E:_ 1 is 


e -iw e iwUH/2) 


In the local description, Uj (f ) = £f =0 Uj L (/). For the calculation of eigenvalues below, we assume the 
data share the above property with the exact solution: u ; _ 1 ( f) = e _/w u / (/), or, for all i, 0 < i < /c, 

uj- lit = e ~ Iw u Jit . (A. 3) 


To calculate the eigenvalues, the /c + 1 values 1 are grouped together as a vector: with superscript T 
denoting the transpose, set 

Uj = (u jj0 ,...,u jjk ) T . (A.4) 

The solution of (A. 1 ) by the methods of the previous sections on each cell Ej can be expressed as 

d 

— uj = + C 0 Uj (A.5) 

where C_ x and C 0 are (/c + 1) X (/c + 1) matrices. Using (A.3), we replace u J _ 1 by e~ IW Uj. The spatial 
discretization results in 


d 

— U, = SUj 

dt ’ ‘ 

where the (k + 1) x (/c + 1) matrix S (for ‘space’ or ‘semi-discrete’) is given by 


(A.6) 


5 = e- /w C_! + C 0 . (A. 7) 

Equation (A.6) is similar to the differential equation du/dt = Au whose solution is ce Xt where c is an 
arbitrary constant. If Re(/L) < 0, the solution is stable. Here, the eigenvalues of 5 take the place of A. For 
all schemes discussed here, 5 has k + 1 eigenvalues. Among these, the principal eigenvalue , denoted by 
S(w), is the one approximating the exact value - Iw of (A.2). All eigenvalues must lie in the left half of 
the complex plane for the semi-discretization to be stable. The collection of all eigenvalues forms the 
spectrum of the scheme. 
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To find the order of accuracy, note that with a uniform mesh of width h, the principal eigenvalue S (w) 
‘approximates’ —hd/dx. A scheme is accurate to order m if S ‘approximates’ —hd/dx to 0(h m+1 ); more 
precisely, for small vv, 

S(w) = —Iw + 0(w m+1 ). (A.8) 

In practice, we obtain the order of accuracy of a scheme by the following procedure. First, set w to be, 
say, 7r/4. We can estimate the error 

E(w) = S(w) + Iw (A.9) 

By halving the wave number w (which is equivalent to doubling the number of mesh points), the error 
corresponding to w/2 (= n/8) is 

E(w/2 ) = S{w / 2) + /w/2. (A.10) 

Since 0((w/2) m+1 ) « (l/2) m+1 0((w) m+1 ), for a scheme to be m - th order accurate, the following 
condition must hold: E(w)/E(w/ 2) ^ 2 m+1 . That is, the order of accuracy is given by 

m ~ L ° g {^k) /Loe(2) - 1 ' <A11) 
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